

/*
 * Class:        GammaProcessPCA
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @authors      Jean-Sebastien Parent and Maxime Dion
 * @since        july 2008

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;
import umontreal.iro.lecuyer.stat.*;
import umontreal.iro.lecuyer.stat.list.*;


/**
 * Represents a <SPAN  CLASS="textit">gamma</SPAN> process sampled using the principal
 * component analysis (PCA).  To simulate the gamma process at times
 * 
 * <SPAN CLASS="MATH"><I>t</I><SUB>0</SUB> &lt; <I>t</I><SUB>1</SUB> &lt; <SUP> ... </SUP> &lt; <I>t</I><SUB>d</SUB></SPAN> by PCA sampling, a Brownian
 * motion 
 * <SPAN CLASS="MATH">{<I>W</I>(<I>t</I>), <I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN> with mean <SPAN CLASS="MATH">0</SPAN> and variance parameter <SPAN CLASS="MATH"><I>&#957;</I></SPAN> is
 * first generated at times 
 * <SPAN CLASS="MATH"><I>t</I><SUB>0</SUB> &lt; <I>t</I><SUB>1</SUB> &lt; <SUP> ... </SUP> &lt; <I>t</I><SUB>d</SUB></SPAN>
 * by PCA sampling (see class {@link BrownianMotionPCA}).
 * The independent increments 
 * <SPAN CLASS="MATH"><I>W</I>(<I>t</I><SUB>j</SUB>) - <I>W</I>(<I>t</I><SUB>j-1</SUB>)</SPAN> of this process
 * are then transformed into independent <SPAN CLASS="MATH"><I>U</I>(0, 1)</SPAN> random variates <SPAN CLASS="MATH"><I>V</I><SUB>j</SUB></SPAN> via
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay">
 * <I>V</I><SUB>j</SUB> = <I>&#934;</I>((&tau;_j-&tau;_j-1)<SUP>1/2</SUP>[<I>W</I>(<I>&#964;</I><SUB>j</SUB>) - <I>W</I>(<I>&#964;</I><SUB>j-1</SUB>)]),&nbsp;&nbsp;&nbsp;&nbsp;<I>j</I> = 1,..., <I>s</I>
 * </DIV><P></P>
 * Finally, the increments of the Gamma process are computed as
 * 
 * <SPAN CLASS="MATH"><I>Y</I>(<I>t</I><SUB>j</SUB>) - <I>Y</I>(<I>t</I><SUB>j-1</SUB>) = <I>G</I><SUP>-1</SUP>(<I>V</I><SUB>j</SUB>)</SPAN>, where <SPAN CLASS="MATH"><I>G</I></SPAN> is the gamma distribution
 *  function.
 * 
 */
public class GammaProcessPCA extends GammaProcess  {
    double[] arrayTime;
    BrownianMotionPCA BMPCA;




   /**
    * Constructs a new <TT>GammaProcessPCA</TT> with parameters
    * 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#957;</I> = <texttt>nu</texttt></SPAN> and initial value 
    * <SPAN CLASS="MATH"><I>S</I>(<I>t</I><SUB>0</SUB>) = <texttt>s0</texttt></SPAN>.
    * The random variables are created using  <TT>stream</TT>.
    * Note that the same {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    *  is used for the <TT>GammaProcessPCA</TT> and for the
    * {@link BrownianMotionPCA} included in this class.  Both the
    * {@link GammaProcessPCA} and the {@link BrownianMotionPCA} are generated by
    * inversion.
    * 
    */
   public GammaProcessPCA (double s0, double mu, double nu,
                           RandomStream stream) {
        super (s0, mu, nu,  new GammaGen (stream, new GammaDist (1.0)));
        this.BMPCA = new BrownianMotionPCA(0.0, 0.0, Math.sqrt(nu), stream);
    }


   /**
    * Constructs a new <TT>GammaProcessPCA</TT> with parameters
    * 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#957;</I> = <texttt>nu</texttt></SPAN> and initial value 
    * <SPAN CLASS="MATH"><I>S</I>(<I>t</I><SUB>0</SUB>) = <texttt>s0</texttt></SPAN>.
    * All the random variables, i.e. the gamma ones and the normal ones,
    * are created using the {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    * included in the {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen}
    * <TT>Ggen</TT>. Note that the parameters of the
    * {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen} object
    * are not important since the implementation forces the generator
    * to use the correct parameters (as defined above).
    * 
    */
   public GammaProcessPCA (double s0, double mu, double nu, GammaGen Ggen) {
        super(s0, mu, nu, Ggen);
        this.BMPCA = new BrownianMotionPCA(0.0, 0.0, Math.sqrt(nu), Ggen.getStream());
    }


   public double[] generatePath() {
        double[] uniformsV = new double[d];
        arrayTime = BMPCA.getObservationTimes();
        int i;
        double[] BMpath = BMPCA.generatePath();
        double sigma;
        for(i = 0; i < d; i++){
            sigma = BMPCA.getSigma() * Math.sqrt(arrayTime[i + 1] - arrayTime[i]);
            uniformsV[i] = NormalDist.cdf01( ( BMpath[i+1] - BMpath[i] )/sigma);
        }
        path[0] = x0;
        for(i = 0; i < d; i++){
            path[i+1] = path[i] +
               GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniformsV[i]);
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   public double[] generatePath (double[] uniform01)  {
        double[] uniformsV = new double[d];
        arrayTime = BMPCA.getObservationTimes();
        int i;
        double[] BMpath = BMPCA.generatePath(uniform01);
        double sigma;
        for(i = 0; i < d; i++){
            sigma = BMPCA.getSigma() * Math.sqrt(arrayTime[i + 1] - arrayTime[i]);
            uniformsV[i] = NormalDist.cdf01( ( BMpath[i+1] - BMpath[i] )/sigma);
        }
        path[0] = x0;
        for(i = 0; i < d; i++){
            path[i+1] = path[i] +
               GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniformsV[i]);
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * This method is not implemented in this class since the path
    * cannot be generated sequentially.
    * 
    */
   public double nextObservation() {
       throw new UnsupportedOperationException ("nextObservation is not implemented in GammaProcessPCA");
    }


   /**
    * This method is not implemented in this class since the path
    * cannot be generated sequentially.
    * 
    */
   public double nextObservation (double nextT)  {
       throw new UnsupportedOperationException ("nextObservation is not implemented in GammaProcessPCA");
    }


   /**
    * Returns the {@link BrownianMotionPCA} that is included in the
    * {@link GammaProcessPCA} object.
    * 
    */
   public BrownianMotionPCA getBMPCA()  {
        return BMPCA;
    }



   /**
    * Sets the observation times of the {@link GammaProcessPCA} and the
    * {@link BrownianMotionPCA}.
    * 
    */
   public void setObservationTimes (double[] t, int d) {
        super.setObservationTimes(t, d);
        BMPCA.setObservationTimes(t, d);
    }


   /**
    * Sets the parameters <TT>s0</TT>, <SPAN CLASS="MATH"><I>&#956;</I></SPAN> and <SPAN CLASS="MATH"><I>&#957;</I></SPAN> to new values, and sets
    * the variance parameters of the {@link BrownianMotionPCA} to <SPAN CLASS="MATH"><I>&#957;</I></SPAN>.
    * 
    */
   public void setParams (double s0, double mu, double nu) {
        super.setParams(s0, mu, nu);
        BMPCA.setParams(0.0, 0.0, Math.sqrt(nu));
    }


   /**
    * Resets the {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    * of the gamma generator and the
    * {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} of
    * the inner {@link BrownianMotionPCA} to
    * <TT>stream</TT>.
    * 
    */
   public void setStream (RandomStream stream)  {
        super.setStream(stream);
        this.BMPCA.setStream(stream);
}

}

